Exterior complex scaling as a perfect absorber in time-dependent problems 



O 
(N 

(N 



Oh' 

I ■ 

Oh- 

a 

o 

q 

O 

Oh. 



> 

o 
in 

(N 

o 
o 



Armin Scrinz{3 

Ludwig Maximilians University, Theresienstrasse 37, 
80333 Munich, Germany, EU 
and 

Wolfgang Pauli Institute, 1090 Vienna 
(Dated: February 12, 2010) 

It is shown that exterior complex scaling provides for complete absorption of outgoing flux in 
numerical solutions of the time-dependent Schrodinger equation with strong infrared fields. This 
is demonstrated by computing high harmonic spectra and wave-function overlaps with the exact 
solution for a one-dimensional model system and by three-dimensional calculations for the H atom 
and a Ne atom model. We lay out the key ingredients for correct implementation and identify 
criteria for efficient discretization. 
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INTRODUCTION 

The absorption of outgoing parts of the wave function 
at the boundaries of a finite volume is a key problem for 
any efficient numerical solution of the time-dependent 
Schrodinger equation (TDSE) and it has been amply 
dealt with also in recent literature (see, e.g., Q and ref- 
erences therein). This interest has been renewed in the 
context of intense laser-matter interactions: speaking in 
terms of physics, strong fields lead to large ionization 
and therefore large fiuxes out of a central region. For 
strong field induced electronic and nuclear dynamics in 
atoms and molecules and high harmonic generation, elec- 
trons far from the system play no role and can be dis- 
regarded. When solving the TDSE for these processes, 
one can therefore identify an inner region (a finite vol- 
ume) where an exact solution is of interest. Outside that 
region one must, by some means, truncate the solution 
without compromising the inner region. This is particu- 
larly important for higher-dimensional problems involv- 
ing two or more electrons in order to control the size 
of the discretization. Out of the large number of ap- 
proaches towards that goal the majority of computations 
of strong laser-matter interactions employed one of the 
following methods: absorbing masks 2j|, complex absorb- 
ing potentials (CAPs) 3, 

and exterior complex scaling 

(ECS) id. 

The two recent numerical studies on ECS have cast 
doubt on the efficiency |i5] and maybe even the funda- 
mental correctness the method in numerical practice 
In the present paper we will show that ECS indeed is a 
perfect absorber to full computational accuracy (14 dig- 
its). In addition, it allows highly efficient implementa- 
tion where only a small fraction of the total discretiza- 
tion points are used for absorption. In both respects it 
far outperforms commonly used monomial CAPs. As a 
third point, as noted early on &\, ECS is not just an ab- 
sorber: ideally, it keeps a record of the dynamics in the 
outer region, which, in principle, could be recovered. We 



will provide numerical evidence for this fact. 

After giving a brief review of ECS, we will present with 
some care our discretization method, as it plays an im- 
portant role for correct and efficient implementation of 
ECS. The general characterization of ECS and a compar- 
ison with CAPs is done using a one-dimensional model 
system, and finally we will present results in three dimen- 
sions for the hydrogen atom and a single-electron model 
of Ne. 



TDSE WITH A LASER FIELD 



We want to solve the TDSE of the general form 



^-*(x,t) = 



^{x,t), (1) 



where x will be either a single x or three x,y,z spatial 
coordinates. Aj', then denote ^r,^ and Laplace 
and Nabla operators, respectively. V{x) is a system- 
dependent binding potential and Alt) is the vector po- 
tential of the laser field. Here we have chosen the ve- 
locity gauge and removed the term Ait)"^ /2 by a time- 
dependent unitary transform. As the initial state we use 
the lowest energy eigenfunction of the field-free Hamilto- 
nian operator — 5A -I- V. We use vector potentials with 
finite duration 



A{t) = ^0 cos^ 
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in the time interval [—nT,nT] with n = 1,5,10. The 
peak vector potential is Aq — Aq and Aq = {Aq, 0, 0) in 1 
and 3 dimensions, respectively. Such pulses with a single 
or a few oscillations of the electric field, linear polariza- 
tion, and peak field amplitude at i = are frequently 
used as models in numerical studies. 

The complete information of the system inside some 
inner region \x\ < Ro is contained in the wave function 
amplitude. For characterizing the accuracy of our results 
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by a single number, we use the overlap between an "ex- 
act" solution Vl/ex obtained from a calculation in a very 
large box and the approximate solution 5* 



1 



l*cx|||| 



(3) 



The scalar product is restricted to the inner region or a 
sub-region of the inner region i? C {|x| < i?o} 



(*cx|*>S 



(4) 



and II • \ \b is the corresponding >C^-norm. 

A quantity of immediate physical interest is the in- 
tensity spectrum of the harmonic response given by the 
Fourier transform of the "acceleration of the dipole" 



5*H= J- 



(*(f,t)||^|*(f,t))s 



n > 2 



(5) 



For the comparison, integrals are restricted to the inner 
region B. In general S{uj) is a highly oscillatory quantity 
varying by several orders of magnitude. The local error 
of the spectrum relative to an "exact" spectrum is 



6uJ[S^lf{uJ) — S^^/^^{UJ)] 
fuj-\-Suj 1 c r \ 



(6) 



Local averaging in the denominator suppresses spurious 
spikes due to near-zeros of the spectrum. 



OUTLINE OF ECS THEORY 

There is a large volume of literature available on com- 
plex scaling in general (see, e.g., 0-0]) and on exterior 
complex scaling in particular (see, e.g., U l3l and refer- 
ences therein) . We restrict our summary to communicat- 
ing the basic idea and to emphasizing the points that are 
essential for correct numerical implementation. For this 
we closely follow earlier work, Ref. ^13.] . In one dimen- 
sion, exterior complex scaling consists in continuing the 
coordinate outside a "scaling radius" Rq into the complex 
plane 



, ( X for |a;| < Rq , , 

'^^^«^"(^) = \e»«(x±i?o)Ti?o for Tx > Ro ' 

The effect of the transformation on plane waves at values 

X > Ro is 

For positive p — outgoing waves to the right side — the 
functions become exponentially damped with increasing 
X, while for negative p they grow exponentially. The 
corresponding situation with reversed signs arises for 



X < —Rq. By complex scaling we can distinguish in- from 
outgoing waves simply by their normalizibility without 
the need to analyze the asymptotic phase. In a typical 
discretization we implicitly or explictly use only square- 
integrable functions, by which we exclude ingoing waves 
from a complex scaled calculation. This is the key to 
complex scaling as a perfect absorber: in a well-defined 
region we have a simple instrument to systematically sup- 
press ingoing waves by just requiring that our solution re- 
main square-integrable. A more elaborate version of this 
reasoning can be found, e.g., in the appendix of Ref. [14 1. 

The mathematically rigorous theory of complex scaling 
often uses the alternative point of view that not the wave 
functions, but rather the operator itself is scaled, while it 
remains an operator acting on an ordinary Hilbert space 
of square-integrable functions. We follow this line of rea- 
soning for pointing out a fact of immediate computational 
relevance. We only give here plausibility arguments and 
refer the reader to Ref. [l^ for a more extensive discus- 
sion and references to mathematical literature. 

One starts from real scaling, i.e. replacing i9 in Eq. ^ 
with a real number A and observes that the transforma- 
tion 



*(j;) 



for x<Ro 



g\/2^(^g\(^xTRo)TRo) for T^>i?o 



(9) 

is unitary. The scaling factor e'*'/^ is essential to ensure 
unitarity. Formally, this transform can just as well be 
applied to the Hamiltonian by defining 



\Ro 



(10) 



It is important to note that if H is defined on differen- 
tiable functions '5, the transformed operator is defined 
on the discontinuous functions '^xrq = Uxro"^ and its 
action on these functions is given by 

HxRo^^XRo - UxRaHU^j^^-^xRo = UxR„H^. (11) 

As a unitary transform Uxro leaves the operator's spec- 
trum unchanged and the scaled dynamics is in a one- 
to-one relation to the unsealed. Now the hard part of 
mathematical theory sets in: for a certain class of "di- 
lation analytic" potentials, the operators Hxrq can be 
analytically continued to complex values X ^ i9 without 
changes in the bound state spectrum ^9'] . The continuous 
spectrum is rotated around the continuum threshold into 
the lower complex energy plane by the angle 26. This is 
trivial to see for the free particle and the case Rq — 
where the spectrum cr(— A) transforms as 



cr(-A) = [0,oo) ^ cr(- 



-2i0 



A) = [0,e-'*''oo). (12) 



This property of the continuous spectrum persists when 
dilation analytic potentials are added and for Rq > 0: the 
complex scaled Hamiltonian retains a distinct "memory" 
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of the unsealed Hamiltonian. Proof of dilation analytie- 
ity ean be difficult to find. Beyond some large Rq, where 
most physical potentials have simple decaying tails, the 
expected ECS properties are confirmed by numerical ex- 
periment. 

One can now write an exterior complex scaled TDSE 



«'ejio(x,t).(13) 



Here it is assumed that the potential can be analytically 
continued to complex values Veiig{x) = V[z0iig{x)]. One 
may hope that the dynamics generated by (jl3p is related 
to the original dynamics, and that for |x| < Rq the solu- 
tion is identical to the unsealed solution "^eRo {x) = vE'(a;). 
We will demonstrate below that this expectation can be 
confirmed by numerical experiment to the level of full 
machine precision. 

The main purpose of this brief discussion of ECS the- 
ory is to stress the importance of the correct discontinuity 
in the wave function for defining differential operators. 
The discontinuity at Rq is intimately linked to the uni- 
tarity of the real scaled problem, which in turn secures 
the conservation of spectral properties of the scaled op- 
erator. For given Rq and it has the explicit form 



^8RoiRo-0) = e''/HgH„{RQ + 0) 



(14) 
(15) 



The discontinuity condition p5|) on the derivative arises 
from transforming the continuous first derivatives of 
the original functions. On such functions, one can 
define the complex scaled Laplacian in analogy to 
Eq. (ITT]) by "back-scaling" the scaled function jfg (x) 
'^BRoi^^^^ix T Ro) ± -Rq)) applying the ordinary Lapla- 
cian, and forward-scaling the result: 



(A(,fl„*e%)(a;) = 



A^^eR^x) 



■0Ro{X) for |a;|<i?o 

e-^^'A^gR„ix) for |.T|>i?o. ^ ' 



The factor e~^*^ appears, as the derivative d^/dx^ is ap- 
plied to the back-scaled function rather than to 'i>gfig{x). 
On continuous functions the scaled Laplacian (and any 
derivative) is not defined as an operator in the Hilbert 
space, just as an ordinary Laplacian is not defined on 
discontinuous functions. 

Finally we want to point to the fact that the adjoint 
operator (Ajo.ej)^ = A[o._e] requires functions with the 
complex conjugate condition of Eqs. (|14I15I) . This means 
that for our discretization by a basis set the disconti- 
nuities must not be conjugated when going from ket- to 
bra-vectors. We will show below how this can be easily 
implemented in a finite element basis. 



DISCRETIZATION 

For the discretization of ECS two points are impor- 
tant: (i) the correct implementation of the discontinuity 
and (ii) good approximation of analyticity. Both can be 
most conveniently accommodated in a finite element dis- 
cretization of high rank. 

We follow the implementation strategy laid out in 
Ref. [l3| : each coordinate axis is divided into N elements 
[xn^i, Xn], n = 1, . . . , N. On each element we choose a 
set of Pn linearly independent functions that can be trans- 
formed to obey the conditions 



A"\xn-l) 



fpn\xn) 



else 



(17) 



We will call Pn the "rank" of the finite element. In prin- 
ciple any set of functions that obeys ([17]) can be used in 
a finite-element scheme. In practice we use real-valued 
polynomials which for enhancing numerical stability we 
transform to 

fi^\x) f'f'\x)dx = TO-"^(5y Vij ^ lpn,Pnl 

(18) 

with normalization constants rnl"-* . For the element func- 
tions (|17p Dirichlet boundary conditions are implemented 
by omitting the functions /}^'' and /p^''. Alternatively, 
on the leftmost and rightmost intervals we use polyno- 
mials times an exponential e^°"^ with -I- and — signs on 
the intervals (— oo,a;i] and [xn-i,oo), respectively. The 
conditions on the end element functions are 



f^'\xi)^0 except £\xi) = I 
fi'^\xN-i)^0 except f[^\xN- 



= 1. 



(19) 



The exponent a can be adjusted for best performance, 
but its exact value was found to be uncritical for ECS. 

The finite-element ansatz for the total wave function 
is as usual 

N p„ 

*(^'0 = EE^i"^w/i"^(^)- (20) 

n— 1 i—1 

By construction of the fj;^^\ Eqs. (|17l) . continuity across 
element boundaries is assured by demanding 



,N. 



n = 2, 

Elementwise overlap and Hamiltonian matrices are 



" f]''\x)dx 



(21) 



J Xn-i 









(22) 

H{i)&\x)dx, (23) 



where J{x) denotes the Jacobian function for integration 
over X. The elementwise matrices are added into the 
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H 



(n-l) 



(n-l) 



(") 



(n) 







(n) 
2p„ 



rr("+l) 
^21 



FIG. 1: Placement of the elementwise block hI"^ in the over- 
all Hamiltonian matrix H 



overall discretized matrices H and S such that the last 
row and column of each elementwise matrix overlaps with 
the first row and column of the following element (see 
Fig- [J), which is equivalent to setting the corresponding 
coefficients equal, Eq. ([2T|). As always in finite element 
methods, continuity of the first derivative does not need 
to be imposed (see Ref. [S] for a more detailed discus- 
sion). H and S are M x Af matrices with 



j:n=lPn-N -Ifov all 



M 



< oo 

Y^n=lPn - + 1 for \xq\ = \xn\ = 



(24) 



For ECS we choose the scaling radii to coincide with 
the element boundaries Xn± = ±i?o- The scaled element- 
wise Hamiltonian matrices are evaluated by substituting 
in ((23)) the Jacobian J{x) and the operator H{t) with 
their ECS equivalents 



in) 

0Ro,i3 



r 



^10 



dxJfl'''>Hff 
(»), 



An) 



\x\ < Ro 
\x\ > Ro 



(25) 

As we use real functions /^^"^ we could omit the complex 
conjugation and the resulting matrices are complex sym- 
metric, i.e Hjj^^^ .j ^ Hg"^^ -^. The discontinuity ([14]) is 
brought into the system by the factor e*^ for the inte- 
grals I a; I > Rq: it amounts to multiplying all functions 
/j^"^ outside the scaling radius by e'^^^ and as the dis- 
continuity does not get complex conjugated, the bra and 
ket discontinuity factors do not cancel but multiply to e'^ . 
Like in the unsealed case, the continuity condition on the 
first derivative does not need to be imposed for fi- 
nite elements. The procedure for constructing the overall 
matrix HgR^ is identical to the unsealed case. Replacing 
H{t) by 1 results in the correct (non-hermitian) overlap 
matrix Sgng for the discretized problem. In practice, the 
matrices Hbr^ and SgRg are rarely set up explicitly, as 



applying the elementwise matrices to the corresponding 
sections of the coefficient vectors is far more efficient. 

There are no specific issues for time propagating the 
discretized system 



(26) 



except maybe that very high accuracy was needed for our 
comparisons. If anything, ECS mitigates the well-known 
stiffness problems for explicit time-integrators, as high 
kinetic energies are also associated with large imaginary 
parts and decay rapidly. We use Runge-Kutta schemes 
with self-adaptive step size and self-adaptive order up 
to order 7. Robust error control is achieved by single-to- 
double-step comparisons. We obtain significant speedups 
of the propagation by removing states with very high 
eigenvalues of the field-free Hamiltonian from the simu- 
lation space by explicit projection. 



ECS FOR A 1-D PROBLEM 

We first investigate ECS for the 1-dimensional "hydro- 
gen atom" with the model potential 



V{x) 



1 



V2- 



(27) 



which gives the ground state energy —1/2. Here and 
below use the peak vector potential \Ao\ — 1.26 and 
the optical period T — 104.8. If interpreted as atomic 
units, these parameters correspond to peak intensity 
2 X lO^^W/cm^ and wave length 760 nm. We will show 
results for FWHM of amplitude of n = 1, 5 and 10 optical 
cycles and total pulse durations of 2,10, and 20 optical 
cycles, see Eq. The classical quiver amplitude of an 
electron in this field is x T/27r = 21 atomic units. At 
the end of a single cycle pulse with this intensity around 
20% and after a 5 cycle pulse more than 80% of the elec- 
tron probability have left the range [-40,40]. 

Within the framework of this model system we will an- 
swer the following questions: Can ECS be considered a 
perfect absorber? Can the scaling radius be put inside 
the range of the quiver motion, i.e. Rq < 21? Does ECS 
work for long pulses? Which parameters determine the 
efficiency of ECS? How many discretization coefficients 
are needed? How does ECS perform compared to mono- 
mial CAPs? Does ECS work for length gauge? 



ECS is a perfect absorber 

We call an absorber perfect, if the error £[—Ro, Rq] de- 
fined in Eq. (jS]) is on the level of machine precision. As 
the "exact" result ^fgx for comparison we use a unsealed 
calculation on a large box [x2,xn-i\ = [—1180,1180] 
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with a total of M = 4801 discretization coefficients dis- 
tributed over 120 elements with constant rank p„ = 41. 
The elements are equidistant except for the infinite end 
elements xq = — oo and xpf = oo with exponent a — 0.5, 

Eq. dm). 

For ECS we use the parameters 9 — 0.5 and Rq — 40 
and finite elements that up to Rq are the same as in the 
unsealed calculation. In the scaled ranges on either end of 
the axis we use infinite elements (cx),— i?o] and [Ro,oo) 
with pi = pn = 41 and exponent a = 0.5. At this 
point, no attempt was made to minimize the number of 
coefficients used for absorption by optimizing the scaling 
parameters. Indeed, with the given parameters we obtain 
for the £^ errors at the end of the pulses t — nT 

F\ R R 1-/2x10"'' forn = l 



The error of the wave function amplitude is about the 
square root of these values and it remains constant after 
the initial rise, see Fig. ([2]). The error level is constant 
through the whole range [~R(),Ro] and there is a sharp 
edge to the scaled region, where the wave function is not 
directly related to the unsealed one. The errors indicate 
the accuracy limits of our numerical integration scheme 
and are not determined by ECS. It is therefore fair to 
say that, at least for the present model, ECS acts as a 
perfect absorber. 




FIG. 2: (color online) Evolution of the relative error |^(a;) — 
^ex(2;)|/i*I'(2:)| during a 5-cycle pulse. The denominator is 
averaged over 5 grid points to avoid spurious spikes. For the 
pulse parameters and discretization, see text. The sharp rise 
of error marks the boundaries of the inner region. A plane is 
drawn at error level 10~^ (color blue); only a few error peaks 
in the inner region are above 10~®. Away from the center, 
relative errors are enhanced initially as the wave function is 
nearly zero. 



Element rank and infinite end elements 

The choice of conspicuously high element rank for 
these very accurate calculations is not by coincidence. 
Complex scaling depends on analyticity properties of the 
Hamiltonian. It is therefor not surprising that we ob- 
serve a strong dependency of the accuracy on the de- 
gree to which our discretization can approximate ana- 
lytic functions. Any localized basis, such as finite ele- 
ments or B-splines is not analytic by definition because 
of the finite support of the basis functions. However with 
increasing polynomial degree, loosely speaking, one gets 
closer to analytic functions. Table U lists the error of the 
wave function in the range [-35,35] for increasing element 
rank. As ±i?o must fall onto element boundaries, we had 
to choose slightly different values Rq for the different el- 
ement ranks. The ECS absorption range is discretized 
with between 36 and 45 exponentially damped functions 
with a = 0.4 such that the sum of coefficients in the 
scaled and unsealed regions was M = 241 for all calcu- 
lations. For the error estimates at each pn a large box 
real calculations was performed with same Pn and the 
same number of points as for ECS in < Rq. From 
Table H] we see that, depending on the desired accuracy, 
it is advisable to use polynomial degrees 8 or higher. 

TABLE I: Dependence of the final wave function error on the 
element rank pn- AH ECS calculations are for a single cycle 
pulse and a total of M = 241 discrete coefficients. 
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For practical purposes we want to mention that the 
variation of an ECS calculation with 9 and box size is 
not a safe indicator of its accuracy. We found ECS cal- 
culations with fixed Rq, element rank and element sizes 
to be far better consistent among each other than their 
error relative to the unsealed result. For reliable accuracy 
estimates one must vary the scaling radius Rq. 

The use of infinite elements at the ends of the simu- 
lation box greatly contributes to the good performance 
of ECS. Table [III compares a few finite-box calculations 
with a calculation using infinite end-elements with only 
21 discretization points. Only at rather large finite boxes 
and a larger number discretization points one reaches the 
infinite element result. 

The explanation for this may be as follows: it was 
noticed in Ref. [HI that long wave lengths cannot be ac- 
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TABLE II: Error of ECS calculations with infinite and fi- 
nite absorption ranges. In all calculations we use a single- 
cycle pulse, rank p„ = 21 and 160 discretization points in 
[—Ro,R(i] ~ [—40,40]. The length of the absorption range is 
^ = iio — a^i = — Ro and A/a is the number of coefficients 
for absorption at each side. 

A Ma £[-40,40] 

oo 21 2 X 10"^^ 

20 20 4 X 10"* 

40 40 2 X 10"" 

60 60 1 X lO"'' 

80 80 1 X 10"^^ 

commodated in a finite ECS region and deteriorate ab- 
sorption by reflections. Such long wave lengths have very 
little structure and should be easily parameterizable. It 
seems that the exponential tail of our end-element func- 
tions is sufficient to accommodate slowly varying long 
wave-length parts of the ECS wave function. We leave 
a more detailed investigation of this observation to later 
work, but conclude that efficient ECS is best done with 
infinite end intervals. 

Choice of Rq and back-scaling 

We find that the quality of the wavefunction in the 
unsealed region is not affected by the choice of the ECS 
radius Rq. Table Hill shows the errors £[—Ro,Ro] for 
Rq = 5, 10, 20 and 40. The general error level in these 
calculations is slightly higher as we used a lower element 
rank of p„ = 11 in order to be able to make the two ele- 
ments of the inner region small. The density of discretiza- 
tion points was kept constant through all calculations. 
We see that the error level is independent of whether 
the ECS radius is chosen inside Rq = 5, 10, 20 or outside 
i?o — 40 the classical quiver amplitude of ao ~ 21. Er- 
rors only start to rise, when the total size of the box indi- 
cated by the number of discretization points M becomes 
too small. This may not be surprising, if we assume that 
the spatial range of the dynamics remains essentially un- 
changed by complex scaling: if the box, be it scaled or 
unsealed, cannot let a particle go the full distance and 
then return without reflections, distortions must occur. 

TABLE III: Dependence of the final wave function error on 
ECS radius Ro and on the total number of discretization 
points M. 

M Ro £[-Ro,Ro] 
241 40. 1.0 X 10"" 
201 20. 5.6 X 10"^^ 
160 10. 2.9 X 10"^^ 
160 5. 1.5 X 10"^^ 
100 10. 1.8 X 10"^^ 
80 10. 1.2 X 10"® 
60 5. 3.6 X 10"^ 



There is an interesting conclusion that one may draw 
from this independence of Rq: the fact that significant 
flux moves into the scaled region and then back out with- 
out corrupting the unsealed part of the wave function 
indicates that also in the scaled region the TDSE dy- 
namics is encoded correctly, although in a different way. 
Our numerical results are a striking corroboration of this 
conjecture that was made early on in ECS theory 
In principle one may hope to recover the unsealed wave 
function by analytic continuation. This hope for back- 
scaling, in fact, was the original motivation for introduc- 
ing the analytic form of functions on the end-elements, 
as an ordinary finite element function cannot be unam- 
biguously analytically continued. We have not further 
pursued this possibility for two reasons: first, with larger 
Pn = Pi and larger scaling angles 9 we encountered se- 
vere numerical problems, as the back-scaled basis func- 
tions become highly oscillatory and cancellation errors 
destroy the reconstruction of the unsealed wave function; 
the second reason is the striking success of ECS with just 
a few points needed for absorption. It is safer and sim- 
pler to just discard the small absorption range and use 
the inner region directly for the evaluation of physical 
quantities. Yet, if for one reason or another, one wishes 
to back-scale a time-dependent ECS wave function, our 
results indicate that such a procedure can be successful. 
One may in that case use a representation of the scaled 
region that is less plagued by numerical problems than 
our exponential basis. 

Choice of scahng angle 6 and exponent a 

Although with sufficiently large absorption range one 
can always obtain perfect absorption independent of scal- 
ing angle 9 and damping exponent a, optimizing these 
parameters in a given situation allows to obtain good 
absorption with very few absorption points. Figure [3] 
shows the error £[—40,40] for n=l and n=5 cycle cal- 
culations with Ma = 10 and 20 absorption points on 
either end of the interval. The exact choice of the pa- 
rameters is not critical for the ~ 20 calculations, 
where full accuracy is reached for rather large parameter 
ranges. As to be expected, the 5-cycle calculation with 
Ma — 10 is most sensitive to 9 and a, but still in a range 
oi9 — 9q±0.1 and a = ao±0.1 around the optimal values 
ao, ^0 ~ 0.3, 0.6 accuracy deteriorates only by 2 orders of 
magnitude to the still acceptable value of 10"*. There is 
a clear anti-correlation between 9 and a, which may be 
explained looking at the oscillatory behavior of the back- 
scaled exponentiallm exp [— are"*^] = sin [a sin 0r] . We 
conjecture that the effective back-scaled wave number 
"f — a sin 9 is the relevant parameter for efficient absorp- 
tion. Correlation between the parameter 7 and 9 nearly 
vanishes and optimization can safely be performed for 
each parameter independently. 
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FIG. 3: (color online) Error f [— 40, 40] as a function of scaling 
angle 6 and exponent ct for n=l and n=5-cycle pulses. For 
the 5 cycle pulse, longer wave length reach the boundaries; 
optimal exponent and scaling angle are smaller and a longer 
(20 point) absorption range significantly increases accuracy. 



Comparison to complex absorbing potentials 

A popular and comparatively straight forward way of 
absorbing outgoing flux are complex absorbing potentials 
(CAPs). The basic idea is to add at the end of the sim- 
ulation box a potential with a negative imaginary part, 
which leads to exponential damping of the wave function. 
In this simplest form, the method can be considered a 
differential form of absorption by mask functions, where 
at preset intervals a certain part of the wave function is 
removed. A fundamental limitation of these methods is 
that they — even in principle — cannot be strictly reflec- 
tionless. The attempt to obtain minimal reflections has 
lead to range of models, partially including real parts into 
the potential and adjusting to specific physical situations 
(see, e.g. 



16|). 



It is beyond the scope of the present work to perform a 
comprehensive study of CAP for the present type of prob- 
lems. Rather, we use the simple and well-investigated 
monomial CAPs Q 



TABLE IV: Accuracy of ECS and CAP for different absorp- 
tion ranges A and number of absorption coefficients Ma- Scal- 
ing angle 9 and absorption strength a for ECS and CAP, re- 
spectively, were optimized. The errors are calculated at the 
end of a single-cycle pulse. 



Method Ma 


A 


6 or (J 


q 


£[-Ro,Ro 


ECS 


21 


oo 


0.6 




2 X lO"^'^ 


ECS 


20 


10 


0.6 




2 X IQ-" 


ECS 


40 


20 


0.5 




1 X IQ-^ 


CAP 


20 


10 


10-^ 


4 


3 X IQ-^ 


CAP 


20 


10 


2 X 10"" 


6 


4 X IQ-^ 


CAP 


40 


20 4 X lO"** 


4 


3 X IQ-* 


CAP 


60 


30 6 X lO"'^ 


4 


1 X IQ-^ 



W{x) 



(29) 



High harmonic spectra 

Although the error f is a meaningful measure for wave 
function accuracy, it cannot be directly related to the er- 
ror of a given observable. Figure 2] shows the accuracy 
of ECS high harmonic spectra of 1- and 5-cycle pulses 
relative to a real calculation. We find errors on the level 
between 10"* and 10"'' and we could not get much bet- 
ter agreement than this irrespective of discretization and 
scaling parameters. Again this error is related to the 
numerical limits of our discretization and propagation 
schemes: the wave function error is ~ 10"^ and the high 
frequency signal is suppressed by 10~* relative to the fun- 
damental peak, making a relative error of the suppressed 
signal of the order 10"'^ quite plausible. Indeed we find 
similar errors when comparing different, but equally ac- 
curate purely real calculations. More disquieting is the 
^ 1% error at the fundamental frequency, which does 
not appear in large box real calculations. We were not 
able to locate the origin of this error: it persists through 
variations of i?o, specific discretizations, different time- 
discretizations, and also for the 3-d hydrogen calculation 
below (cf. Fig. [5]). The error appears to be related to an 
artificial overall modulation of the signal by the driving 
field, possibly related to internal normalizations during 
propagation. Note that by construction normalization 
errors do not appear in the wave function accuracy mea- 
sure Eq. ([3]). We believe, however, that this error is 
acceptable for all practical purposes. 



for polynomial degrees g = 4, 6 with optimized a in each 
calculation. The criterion for our comparison with ECS 
is the number of discretization points needed for a given 
level of absorption. 

Results are shown in Table HVl With a finite absorp- 
tion range, ECS outperforms CAP roughly by one or two 
orders of magnitude. However, when we use infinite end 
elements with ECS (discretized by only 21 points), we 
can reach absorption to machine precision. We could not 
find a similar adjustment for CAP. 



ECS fails in length gauge 

For field-interaction in length gauge 

dA 

iA{t)-V3^x-— (30) 
at 

ECS completely fails in the time-dependent case. The 
reason for this behavior was pointed out in Ref . Q : when 
length gauge Volkov solutions are complex scaled their 
asymptotic behavior becomes dependent on the sign of 
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TABLE V: Parameters for the Ne model potential Eq. ((SU 



Harmonic order 



FIG. 4: (color online) High harmonic power spectrum for a 
5-cycle pulse (upper panel). Lower panel: Accuracy ^{uj) of 
the high harmonic spectrum with different ECS parameters 
and discretizations. Curve A is the error for {Ro, Ma, 9, ct) ~ 
(40,20,0.7,0.7) relative to a fully converged real calculation. 
The choice of _Ro has has the largest influence on accuracy: 
curve B, the difference between two calculations with _Ro ~ 
40 and Ro ~ 50 closely follows the overall error curve A. 
At fixed -Ro the influence of the other ECS parameters and 
discretization is small: curve C compares a calculations using 
(iZo, Ma, 6*, a) = (40,20,0.7,0.7) and finite element rank p„ = 
21 with (40,40,0.5,0.3) and rank 41. 



the field strength and alternates between damping and 
growth. The convenient distinction between in- and out- 
going waves by their norms is lost. In the language of 
mathematical theory, x is not a dilation analytic po- 
tential, and severely so: complex scaling transforms the 
spectrum of the Stark problem from purely continuous 
into purely discrete and all discrete eigenvalues of the 
scaled Stark Hamiltonian have imaginary parts ■ This 
is sharp contrast to dilation analytic potentials where the 
bound state energies remain unchanged and the contin- 
uous spectrum is only rotated into the lower complex 
plane. 



CALCULATIONS FOR H AND MODEL NE 

In order to demonstrate the apphcability of ECS to 
realistic problems, we show calculations for the H atom 
with 



(31) 



and a single electron model of the Ne atom with the po- 
tential 



1=1 



exp[-Ci|a:| 



(32) 



We use the parameters given in Table |Vl for which our 
model reproduces the ground and first few excited state 



ai Ci 

1 -1 

2 -0.3 0.5 

3 -2.05 2 

4 1.23 1 



energies of Ne. We use linearly polarized pulses with the 
same pulse shape and peak intensity as in the preceding 
section and pulse durations of 1 and 10 optical cycles. 
The calculations are done in polar coordinates with a 
spherical harmonics basis on the angular coordinates and 
high rank finite elements on the r-coordinate. Again an 
infinite last element is used. 

There are no surprises: convergence patterns and ac- 
curacy are very similar to the one-dimensional model. 
Figure [5] shows the harmonic spectrum for hydrogen at 
1 cycle together with errors for different ECS and dis- 
cretization parameters. Error estimate here is by com- 
parison to an Rq — 80 calculation. 




Harmonic order 



FIG. 5: (color online) High harmonic power spectrum from 
the hydrogen atom for a 1-cycle pulse (upper panel). Lower 
panel: Error 'D{iu) with ECS parameters (Ro, AlA,d, a) = 
(40,20,0.5,0.5) relative to a iio = 80 calculation (curve A). 
The relative difference to a calculation with (40,40,0.4,0.4), 
curve B, underestimates the error. The calculation is con- 
verged with 20 angular momenta on the given level of accu- 
racy. More angular momenta do not change the result. At 15 
angular momenta (curve C) accuracy deteriorates. 

No new problems appear due to the more general Ne 
model potential (1321) . FigurelHlshows high harmonic spec- 
tra from a H and Ne for a 10-cycle pulse. Accuracy esti- 
mates were obtained by varying ECS radius Rq. 



DISCUSSION 

As we find high numerical stability and excellent per- 
formance of ECS as an absorber, the questions arises 
what are the reasons for the numerical problems reported 



9 




10 20 30 40 

Harmonic order 



FIG. 6: (color online) High harmonic power spectra from 
the hydrogen and a neon model with a 10-cycle pulse (up- 
per panel). Relative accuracies shown in the lower panel are 
somewhat poorer with the longer pulse, in particular for Ne 
where due to the higher ionization potential the signal is very 
weak. 

in Refs. 0, @, where ECS was applied to essentially the 
same systems. One obvious source of inaccuracies may 
he in possibly low order discretizations. Unfortunately, 
in neither publication an investigation of the dependence 
of the observed effects on discretization is shown. 

Certainly the choice of finite box-sizes lowers the per- 
formance, but according to Table HIl with the very large 
absorption ranges of 80 Bohr used in Ref. Q, excellent 
results should be achievable. 

A possible source of the observed difficulties may be 
the treatment of the overlap matrix. As discussed above 
we replace the ordinary overlap by the pseudo-overlap 
matrix SgRg. With this choice and as we use strictly 
real finite element functions we obtain complex symmet- 
ric matrices {Hgno)'^ = Hgn^ for zero field Aq = Q and 
(SeRo)^ = SsBj,. There are no explicit statements about 
SgRg in Refs. [HQ. Usually, finite difference methods im- 
ply (an approximation to) the identity operator for over- 
lap. The B-spline method used in Q requires a choice for 
SgR^-, and Eq. (20) of Ref. Q seems to imply that indeed 
the identity was used as an overlap matrix. 

The comment on the non-orthogonality of the_eigen- 
vectors of the non-normal scaled Hamiltonian in fs!] also 
seems to indicate, that the ordinary, unsealed overlap 
matrix S was used. Clearly, the eigenvectors fe*^"-' of the 
eigenproblem 

iJeflo^^"^ = Sb^"^E^ (33) 

will not be orthogonal in general. However, we find that 
all eigenvectors of the complex scaled generalized eigen- 
problem 

HgR,C^''> = SgR,(^'^E, (34) 



are pseudo-orthogonal and can be normalized in the sense 

^ — ' \ / Im 

Im 

Then the matrix Hor^-^ has a diagonal representation 

HgR„^Y.^'^^^{^'^T (36) 

i 

and the spectral values Ei appear as discrete approxima- 
tions to the true ECS spectrum with strictly non-positive 
imaginary parts ImEa < 0. We do not have mathemati- 
cal proof for this property of the discrete complex scaled 
system, but we find it vaUd in all our calculations on the 
level of computational accuracy. If on the other hand we 
use the ordinary overlap matrix S, we invariably obtain 
a few eigenvalues Ea with ImEa > which will cause 
long-term instability of the time-propagation. Possibly, 
this is the ultimate reason for the numerical instabilities 
observed in Refs. [3, 

SUMMARY 

We have demonstrated that ECS can serve as a per- 
fect absorber of outgoing flux in the sense that in the 
unsealed inner region it exactly matches a purely real 
calculation on a sufficiently large grid. We were able to 
push the agreement to relative £^ error of lO"""^^. The 
corresponding errors in the wave function amplitude are 
^ 10^^. Both errors are at the limits of our numerical 
integration scheme. 

Furthermore we have evidence that also in numerical 
practice ECS does not just act as an absorber but con- 
serves dynamical information during excursions into the 
absorbing region: even when the quiver motion takes fiux 
deeply into the "absorbing" region, the returning flux is 
identical to the flux in a purely real calculation. 

For this, we found the following points essential: 

(i) implementation of the correct scaled derivatives, in- 
cluding bra-functions with unconjugated disconti- 
nuity, 

(ii) the use of "infinite" absorption ranges, which we 
discretized by polynomials times an exponential, 

(iii) the use of high rank discretization also in the inner 
region to reach the highest accuracies. 

Point (i) leads to a complex symmetric, in particular not 
positive definite discrete overlap matrix which must not 
be approximated by a positive definite matrix. 

Following these rules, we encountered no numerical dif- 
ficulties in the inner region or in the absorbing region, 
using a standard explicit Runge-Kutta scheme for time 
integration. As a tendency, large scaling angles favor 
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good absorption, in many cases we used 9 = 0.7 w 40°, 
which corresponds to an almost purely imaginary contin- 
uous energy spectrum [0, e^'^'^oo). In our basis we found 
the scaling angle ultimately to be limited by numerical 
instabilities due to the complex symmetric overlap ma- 
trix. As excellent absorption can be achieved with as few 
as 20 discretization coefficients in the absorbing region, 
pushing the scaling angle to the numerical limits is not 
necessary in general and scaling angles of^ = 0.3 ^ 0.5 
served well for our purposes. In general, we found the 
scheme numerically robust and not very sensitive to the 
scaling parameters. The option of back-scaling the solu- 
tion to = was abandoned due to severe cancellation 
errors in the related transformations. 

When judging the accuracy of an ECS calculation, it 
is important to vary the ECS radius i?o- Our comparison 
with a real calculation indicates the variation of the result 
with different i?o gives realistic error estimates. Other 
parameters such as rank of the finite elements, length of 
the absorption range, or scaling angle are of secondary 
importance. 

ECS in the present implementation outperforms sim- 
ple monomial CAPs. ECS errors were one or two order 
of magnitude smaller than CAP errors with the same 
spatial discretization. When using infinite end intervals, 
the advantage of ECS can even reach 12 orders of mag- 
nitude! We are aware of the fact CAPs can be greatly 
improved by a variety of measures (see, e.g., (iBl)- How- 
ever, in general these require tuning of the CAP param- 
eters to a given situation. Even with that, we do not 
expect to reach comparable accuracies with CAPs as we 
could demonstrate for ECS. 

We believe that ECS solves the absorption problem 
for the present class of systems in any discretization, 
where implementation of (i)-(iii) is possible. Recovery 
of asymptotic information, such as electron spectra, may 
be inherently difficult as the total amount of informa- 
tion that is contained in the scaled discretization is too 
small, which manifests itself in the ill-conditioning of the 
back-scaling problem. However, in Ref. [Tt! ] we presented 
a scheme for computing electron spectra under the as- 
sumption of a perfect absorber, at that point formulated 



for CAPs. An adaptation of that scheme to ECS and 
extension to few-body dynamics will be investigated in 
future work. 
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